open(BUFF,"/data1/bsi/BORA_processing/devel/eqtl/outfile_plink");
open(BUFF1,"/data1/bsi/BORA_processing/devel/eqtl/outfile_expr_plink");
$line = <BUFF>;
chomp($line);
$line =~ s/^\s+//;
@geno = split(" ",$line);
$line = <BUFF1>;
chomp($line);
$line =~ s/^\s+//;
@expr = split("\t",$line);
for($i=0;$i<@geno;$i++)
{
	$g{$geno[$i]} =1;
}
$k =0;
for($i=0;$i<@expr;$i++)
{
	if(exists($g{$expr[$i]}))
	{
        	$h{$expr[$i]} =$k;
		$k++;
	}
}
for($i=0;$i<@geno;$i++)
{
	if(exists($h{$geno[$i]}))
	{
		$jj[$i] = $h{$geno[$i]}+1;
		$jj2[$h{$geno[$i]}] = $geno[$i];
	}
	else
	{
		$jj[$i] = -1;
	}
}


for($i=0;$i<@expr;$i++)
{
	if(exists($h{$expr[$i]}))
	{
		$jj1[$i] = $h{$expr[$i]}+1;	
		$jj3[$h{$expr[$i]}] = $expr[$i];
	}
	else
	{
		$jj1[$i] =-1;
	}
}
$line_geno = join(" ",@jj2);
$line_expr = join("\t",@jj3);
open(WRB,">/data1/bsi/BORA_processing/devel/eqtl/outfile_result_plink");
print WRB "GENE\tSNP\tA1\tA2\tFRQ\tINFO\tBETA\tSE\tP\n";
#$line = <BUFF>;
#$line = <BUFF>;
#print "@jj\t@jj1\n";
while($line = <BUFF>)
{
	open(WRBUFF1,">/data1/bsi/BORA_processing/devel/eqtl/temp_plink1/expr.fam");
	open(WRBUFF,">/data1/bsi/BORA_processing/devel/eqtl/temp_plink1/expr.dat");
	open(WRBUFF2,">/data1/bsi/BORA_processing/devel/eqtl/temp_plink1/expr.pheno");
	@jaj = split(" ",$line_geno);
	$line_expr = "SNP A1 A2";
	for($y=0;$y<@jaj;$y++)
	{
		$line_expr = $line_expr." $jaj[$y] $jaj[$y]";
		 
	}
	print WRBUFF $line_expr."\n";
	chomp($line);
	@array = split("\t",$line);
	for($i=1;$i<@array;$i++)
	{
		@gg = split(/ /,$array[$i]);
		$gg1 = shift(@gg);
		$gg2 = shift(@gg);
		$gg3 = shift(@gg);
		$gg4 = shift(@gg);
		for($k=0;$k<@gg;$k++)
		{
			if($jj[$k] != -1 )
			{
				$gg1[$jj[$k]-1] = $gg[$k];
			} 
			if($jj[$k] == 0)
			{
				die "$jj[$k] $k\n";
			}
		}
		for($y=0;$y<@gg1;$y++)
		{
			if($gg1[$y] >= 0 && $gg1[$y] <= 0.5)
			{
				$gg1[$y] = 0;
			}	
			elsif($gg1[$y] >= 1.5 && $gg1[$y] <= 2)
			{
				$gg1[$y] = 2;
			}
			else
			{
				$gg1[$y] = 1;
			}
	
		}
		$line = join(" ",@gg1);
		if($gg4 eq '-')
		{
			$gg4 = $gg3;
			$gg4 =~ tr/ATGC/TACG/;
		}
		$line = $gg1." ".$gg3." ".$gg4." $line";
		print WRBUFF $line."\n";
	}
	@array1 = ();
	#$line = <BUFF1>;
	$line = <BUFF1>;
	chomp($line);
	@ee = split("\t",$line);
	$ee1 = shift(@ee);
	for($k=0;$k<@ee;$k++)
        {
        	if($jj1[$k] != -1)
               	{
                	$ee1[$jj1[$k]] = $ee[$k];
                }
        }
	for($y=0;$y<@jaj;$y++)
        {
                print WRBUFF1 "$jaj[$y] $jaj[$y] 0 0 -9 1\n";
		print WRBUFF2 "$jaj[$y] $jaj[$y] $ee1[$y+1]\n";
        }
	system("/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/plink/plink --dosage /data1/bsi/BORA_processing/devel/eqtl/temp_plink1/expr.dat  format=1 --fam /data1/bsi/BORA_processing/devel/eqtl/temp_plink1/expr.fam  --pheno /data1/bsi/BORA_processing/devel/eqtl/temp_plink1/expr.pheno --assoc --out  /data1/bsi/BORA_processing/devel/eqtl/temp_plink1/kk");
	open BUB,"/data1/bsi/BORA_processing/devel/eqtl/temp_plink1/kk.assoc.dosage" or die "gene:$ee1\n";
	$line = <BUB>;
	while($line=<BUB>)
	{
		chomp($line);
		$line =~ s/\s+/\t/g;
		print WRB $ee1."\t$line\n";
	}
	system("rm /data1/bsi/BORA_processing/devel/eqtl/temp_plink1/*");	
}
